El lector podrá ahora darse cuenta por sí mismo que hay diferencias al programar en paralelo, con respecto a la forma usual. La principal de ellas es que el GPU no es un procesador serial, sino un procesador stream. Un procesador serial, también conocido como arquitectura de von Neumann, ejecuta instrucciones serialmente, actualizando la memoria a medida que avanza. Un procesador en stream funciona de una manera ligeramente diferente, ejecutando una función (tal como un fragmento de programa) en un conjunto de registros de entrada, produciendo un conjunto de registros de salida (e.g. pixeles sombreados) en paralelo. Los procesadores en stream, típicamente se refieren a dicha función como kernel y al conjunto de registros como stream (muy imaginativos). La información es lanzada al procesador, operada vía una función kernel, y sacada a la memoria, como se muestra en la figura 1. Cada elemento pasado al procesador es procesado independientemente, sin ser influido por los demás elementos. Esto permite a la arquitectura ejecutar un programa en paralelo sin la necesidad de exponer las unidades paralelas o cualquier construcción en paralelo al programador. Es por dicha razón que al escribir la función kernel en el capítulo pasado no notamos el paralelismo, no fue sino hasta el momento de lanzar el kernel que vimos esto.
Y qué tal si quisiésemos hacer física con esto, sería increíblemente incómodo que el único tipo de dato que pudiésemos manejar fuera el sombreado de un pixel; sin embargo, gracias a los avances en GPU se pueden manejar registros que no necesariamente son pixeles y sombreados; bien podría ser un conjunto arbitrario de datos, por ejemplo puntos de una red y ecuaciones fisicas. Digamos, que queremos implementar una simulación simple de partículas libres y dejar que el GPU lleve a cabo la mayor parte de las operaciones físicas. Usando texturas de punto flotante y buffers de pixel (pbuffers), guardamos la posición de las partículas, valocidades, y orientaciones. Un kernel hace los cálculos necesarios para obtener la nueva posición de cada partícula. Para calcular todas y cada una de las nuevas posiciones simplemente tomamos un cuadrante que tenga tantos pixeles como partículas, tal y como se muestra en la figura 2. Las coordenadas de textura para el cuadrante indican al kernel cuál de las partículas guardadas procesará. El resultado guardado en el pbuffer contiene los valores de las nuevas posiciones.
Este ejemplo de un sistema de partículas demuestra una obvia pero útil aplicación del GPU para el cómputo de propósito general (GPGPU). La operación de actualizar la posición de las partículas puede generalizarse al proceso de aplicar una función a un arreglo. Esta operación - también llamada mapeo en lenguajes de programación funcional - puede ser usada para una amplia variedad de aplicaciones.
Pero ¿qué es este mapeo del que tanto hablamos?
En el mapeo, una función opera en cada elemento de un arreglo de entrada y produce un resultado, por ejemplo sumar a cada entrada de un vector un número, o multiplicar todo el vector por una constante, con la conveniencia de que lo que se lee en la posición i de la entrada sirve para obtener el resultado que se escribe en la posición i de la salida. Para procesar todos los elementos de la entrada en paralelo, una función de mapeo debe ser pura. Esto quiere decir que genera exactamente el mismo resultado para la misma entrada, y que su ejecución no tiene efectos secundarios.
Dado que no hay necesidad de sincronizar a los threads, y no se comparte información, el patrón de mapeo se ajusta perfectamente a arquitecturas paralelas, de muchos núcleos. En las implementaciones en paralelo de los patrones de mapeo, cada thread ejecuta una instancia de una función mapeo y genera su correspondiente resultado. Éste patron es usado en muchos dominios, incluyendo el procesamiento de imágenes, simulaciones financieras, simulaciones físicas, etc.
A continuación mostramos un código para multiplicar un vector de 100 entradas (predefinidas del 1 al 100) por un escalar.
#include "stdio.h"
#define N 100
__global__ void add(int *a, int *c)
{
int tID = blockIdx.x;
if (tID < N)
{
c[tID] = 3*a[tID];
}
}
int main()
{
int a[N], c[N];
int *d_a, *d_c;
cudaMalloc((void **) &d_a, N*sizeof(int));
cudaMalloc((void **) &d_c, N*sizeof(int));
// Llenar el arreglo
for (int i = 0; i < N; i++)
{
a[i] = i;
}
cudaMemcpy(d_a, a, N*sizeof(int), cudaMemcpyHostToDevice);
add<<<N,1>>>(d_a, d_c);
cudaMemcpy(c, d_c, N*sizeof(int), cudaMemcpyDeviceToHost);
for (int i = 0; i < N; i++)
{
printf("3*%d = %d\n", a[i], c[i]);
}
return 0;
}
In [1]:
!nvcc -o ./Programas/multEscalar ./Programas/multEscalar.cu
In [2]:
!./Programas/multEscalar
En este punto la pregunta en el aire debe ser ¿Qué es #define N 100? Bueno, simplemente definimos una constante que se llama N y que tiene como valor 100.
Otra cosa que llama la atención es que en el kernel tenemos esta parte
if (tID < N)
{
c[tID] = 3*a[tID];
}
Por el momento sólo explicaremos que esto se hace porque planeamos lanzar más de 100 threads. "¿Por qué?" es una pregunta que contestaremos más adelante. Por el momento sólo se necesita saber que ese if nos permite hacer trabajar solamente a los 100 primeros threads y deja a los otros sin hacer nada. ¿Qué pasaría si un thread con un índice mayor a 100 intenta ejecutar esto c[tID] = 3*a[tID]?
Visto esto, se deja como ejercicio al lector hacer un código (en la celda de abajo) que sume 3 a cada una de las entradas de un vector (pueden usar el vector que quieran, en especial el que hemos venido utilizando).
In [ ]:
Los patrones dispersar (scatter) y reunir (gather) son similares a los patrones mapeo pero sus accesos a la memoria son aleatorios. Dispersar es una función mapeo que escribe en posiciones aleatorias, y reunir es la combinación de una función mapeo con accesos a memoria que leen de entradas aleatorias. Las implementacionesen paralelo de los patrones dispersar y reunir son similares a las implementaciones del mapeo. Dicho patrón se encuentra comunmente en aplicaciones estadísticas.
Cuando una función combina todos los elementos de un arreglo de entrada para generar una salida única, se dice que está realizando una reducción. Si la función usada por el patrón de reducción es asociativa y conmutativa, por ejemplo, el XOR, el orden en el cuál la operación de reducción es aplicada a las entradas no tiene importancia. En este caso, las implementaciones basadas en árboles pueden ser utilizadas para paralelizar dicha reducción. Las reducciones se pueden encontrar en diferentes áreas, como el aprendizaje de máquinas, física, y estadística; entre otras.
Un ejemplo trivial de esto sería el promediar todas las entradas de un arreglo.
Ejercicio: Escribir un kernel que calcule dicho promedio
In [ ]:
Discuta consigo mismo, o con un amigo, ¿por qué podría no resultar buena idea realizar el promedio en el GPU?
Podemos obtener ayuda de nuestros amigos en los foros de NVIDIA
El escaneo aplica una función asociativa a un arreglo de entrada y genera un arreglo nuevo. Cada n-ésimo elemento del arreglo de salida es el resultado de aplicar la función escaneo a los primeros N (escaneo inclusivo) o N - 1 (escaneo exclusivo) elementos de entrada. El patrón de escaneo es usualmente utilizado en el contexto de procesamiento de señales, y aprendizaje de máquinas.
Discuta un ejemplo en donde se utilice el escaneo
En un patrón "stencil", cada elemento de salida es calculado aplicando la función a su correspondiente elemento en el arreglo de entrada, y a sus vecinos. Éste patrón es útil por ejemplo para hacer borrosa una imagen; si decidimos que en nuestra nueva imagen cada pixel tendrá el color del pixel que estaba originalmente en ese lugar promediado con el color de sus vecinos. Este patrón es usualmente encontrado en procesamiento de imágenes y física.
¿Habrá algún juego divertido en el que el patrón stencil sea útil?
El patrón de partición (o "tile") es similar al patrón stencil. El arreglo de entrada es separado en particiones y cada partición es procesada de manera separada. Cada partición es completamente independiente de las otras. Particionar es utilizando comunmente en aplicaciones que tienen paralelismo de datos para utilizar eficientemente las arquitecturas intrínsecas de jerarquía de memoria y así mejorar el rendimiento. Éste patrón es utilizado comunmente en dominios como procesamiento de imágenes, procesamiento de señales, y modelación en física. Más adelante en las notas se mostrará un ejemplo para exhibir la importancia de este patrón de comunicación.
Esta pregunta no nos llega a la mente en primera instancia, por el hecho de que estamos demasiado acostumbrados a que en la programación serial el orden de ejecución es el orden que se muestra en el código, la computadora no tiene mucho control sobre eso. En programación en paralelo esta cuestión nos preocupa, ¿en qué orden se ejecutarán los threads?
Bueno, veámoslo. El siguiente kernel simplemente imprime el índice del bloque en que está el thread que se ejecuta.
#include <stdio.h>
#define NUM_BLOCKS 16
#define ANCHURA_BLOCK 1
__global__ void hola()
{
printf("Hola mundo! Soy un thread en el bloque %d\n", blockIdx.x);
}
int main(int argc,char **argv)
{
// lanzar el kernel
hola<<<NUM_BLOCKS, ANCHURA_BLOCK>>>();
// forzar a los printf() para que se muestren
cudaDeviceSynchronize();
printf("Eso es todo amigos!\n");
return 0;
}
Lo compilamos, esta vez con la bandera -arch sm_20 para especificar que compilamos en esta arquitectura, esto nos permite compilar sin problemas con la función printf. Se le sugiere al lector intentar la compilacion sin dicha bandera para que se familiarice con este tipo de errores de compilación.
In [3]:
!nvcc -arch sm_20 -o ./Programas/orden ./Programas/orden.cu
In [4]:
!./Programas/orden
In [5]:
!./Programas/orden
In [6]:
!./Programas/orden
Como podemos ver no se ejecutan en un orden predefinido, y ciertamente no lo hacen siempre en el mismo orden. ¿Esto puede ser problemático? Sí, sin embargo más adelante (en el capítulo 5) veremos una forma de lidiar con este problema.